Direct Minimization Generating Electronic States with Proper Occupation Numbers 
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We carry out the direct minimization of the energy functional proposed by Mauri, Galli and Car 
to derive the correct self-consistent ground state with fractional occupation numbers for a system 
degenerating at the Fermi level. As a consequence, this approach enables us to determine the elec- 
tronic structure of metallic systems to a high degree of accuracy without the aid of level broadening 
of the Fermi-distribution function. The efficiency of the method is illustrated by calculating the 
ground-state energy of C2 and Si2 molecules and the W(110) surface to which a tungsten adatom 
is adsorbed. 
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In the Kohn-Sham (KS) formulation of density- 
functional theory , the energy functional of a system 
of electrons is stationary at the point of the ground-state 
energy with respect to the variation of an arbitrary set 
of single-particle orthonormal wave functions {ipi} and 
occupation numbers {rii} where rij varies between zero 
and one. The Janak theorem implies that the occu- 
pation numbers of the ground state at zero temperature 
are chosen such that n, = 1 for the states of KS eigen- 
values below the Fermi level, and = for those above 
the Fermi level. For systems with degeneracies at the 
Fermi level such as metallic systems, however, the deter- 
mination of the occupation numbers of the degenerate 
states is not trivial |^|. One popular approximation 
has been to determine the occupation numbers of the 
states near the Fermi level using the Fermi-distribution 
(FD) or its substitute function where a fictitious temper- 
ature, a parameter to adjust the broadening of levels, is 
included |^||- Further, the thermodynamic free energy 
is minimized instead of the electronic total energy. This 
procedure has the advantage that the calculations of elec- 
tronic wave functions for metallic systems typically con- 
verge more stably and rapidly, and that the broadening 
of each level roughly imitates larger systems or a better 
sampling in the Brillouin zone. In principle, however, no 
physical significance is given to temperature or entropy 
in the context of the zero-temperature KS scheme. It 
is even more serious in practice that theoretical results 
depend on the chosen values of the fictitious broadening 
temperature, and that all of the occupation numbers of 
the degenerate states at the Fermi level are restricted to 
be numerically equal, which is usually incorrect for acci- 
dental degeneracies. 

In this Letter, we demonstrate that the direct min- 
imization (DM) of the energy functional proposed by 
Mauri, Galli and Car (MGC) @ leads to the assign- 
ment of the proper occupation numbers to both degen- 
erate and nondegenerate states; moreover, the DM pro- 
cedure yields the correct self-consistent solutions of the 
KS equation without usage of conventional self-consistent 



field techniques. Consequently, the electronic structure 
of metallic systems at zero temperature can be obtained 
to a high degree of accuracy with no aid of a FD broad- 
ening. We combine the DM method with the real-space 
finite-difference method utilizing the timesaving double- 
grid technique ||. The real-space calculation methods 
have tackled serious drawbacks of the plane-wave ap- 
proach, e.g., its inability to describe strictly nonperiodic 
systems such as clusters and solid surfaces. 

The MGC energy functional for an TV-electron system 
is written as 
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where {</>f } is an arbitrary set of M linearly independent 
overlapping wave functions with spin index <t, which are 
assumed here to be real functions for simplicity, and M is 
taken to be not smaller than the number of the occupied 
states for each spin. F[p\ p^] is the sum of the external, 
Hartree, and exchange-correlation potential energy func- 
tionals, r\ is the electronic chemical potential, Q a is an 
(M x M) matrix: Q\- = 2<5 y - S^, and S a is the overlap 
(</>f |c^J). The charge density is defined as 



matrix: Sf, = 



M 



P(r) = 



<t=U 



CT (r) and p°(r) = £ Q£#(r)tf (r). (2) 



The form of the energy functional Eq. (Q) was originally 
introduced for computation with linear system-size scal- 
ing O(N), and each wave function <f>? was approximated 
to be a Wannicr-likc function localized in an appropri- 
ate region of space (i.e., a localized orbital) to reduce the 
amount of computation. In this paper, we adhere to ordi- 
nary extended wave functions, being free from the errors 
caused by the localization of wave functions. 

E[{4>}, i]] is minimized by a steepest-descent or 
conjugate-gradient algorithm without constraint of the 
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orthogonalization of wave functions. The derivative of 
the functional with respect to the function <f>° is required 
for the minimization, which is given by 



with Tfa 



se[{4>)M 



M 



^[(^s[{4>}]-vM)Q°j 



(3) 



where -H^-s [{</>}] is the KS Hamiltonian. It is straight- 
forward to verify that the ground-state energy is a sta- 
tionary point of E[{ 4>},t)]. Consider the following set 
of the ground-state wave functions {</>f }: 0f°) = 
a i\Xi)j where |xf) is the normalized eigenfunction of 
Hks at {^D = {C } with the eigenvalue ef, i.e., 
HK S [{^f°}]\Xi) = eflxf)) and the coefficient af is set 
such that K| = 1 for ef < r), < \af | < 1 for ef = r), 
and af = for ef > 77. Obviously, {0f } is a stationary 
point of E[{<j)},ri], since 

<5£;[{0},7 7 ]/ ( 5C| Wf}={0f o } =4 a f(l-| a f| 2 )(ef-r,)| x f} 

= 0. (4) 

We now show that the MGC energy functional Eq. (|l|) 
is identical to the functional in the standard form 

M 

E[M,V]= E E<WI-2 V2 ^) 
<r=U • 

+F[p1,p L ]+r l !.N - J p(r)drj, (5) 

with 
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p(v) = E and p CT (r) = |^(r)| 2 , (6) 

<T=T1 

where {ipf } and {nf } are sets of orthonormal wave func- 
tions and the corresponding occupation numbers, respec- 
tively, but nf varies in the range that — 00 < nf < 1. The 
proof is as follows: A set of {<pf } in Eqs. (||) and (||) is 
transformed to an orthonormal set {V'f } using an orthog- 
onalization algorithm, e.g., the Gram-Schmidt method, 
and then ^>f is represented as </>f = Yli ' c Li^i '■ Substi- 



tuting this expansion into Eqs. (0) and (gj, we obtain 
for the energy functional 
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and a similar expression for the charge density, where 

M M 
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cf i . Here, the matrix T CT , of which 
the element of the k row and the I column is , is a non- 
negative definite Hermitian matrix, which can be diago- 
nalized by a unitary matrix U a as (U^T a U <T ) k i = 5 kl \% 
with Af being a non- negative eigenvalue of T a . Hence, 



(9) 



Defining ipf and nf as ip? = Ylf "^j^ji an d n 1 = 
(2 - Af )Af , we obtain Eq. (g) from Eq. (g), and the 
inequality -00 < nf < 1. (Q.E.D.) 

It is very important that the occupation numbers de- 
fined above do not exceed one, because this fact means 
that in the course of the minimization of the the MGC en- 
ergy functional Eq. ([j]) , the Pauli principle automatically 
works to prevent more than one electron from falling into 
each single-particle level. The unphysical negative occu- 
pation gives rise to no severe problem in practical cal- 
culations, as long as initial wave functions for iterations 
are chosen to be close to the ground-state solution. In 
this case, the occupation numbers maintain non-negative 
values during iterations, which is assured by experience 
(e.g., see Fig. |). 

The overall computational scaling in our DM pro- 
cedure combined with the real-space finite-difference 
method of Ref. (|] amounts to 0(M 2 N gri d) operations 
in the calculations of the energy functional Eq. ([[]) and 
its derivative Eq. (g), since the <f>f 's are now assumed 
to be extended wave functions. Here, N gr id is the num- 
ber of coarse grid points in real space. It is noted that 
0(M 2 N gr id) is equal to the scaling order in the orthogo- 
nalization of wave functions indispensable in the conven- 
tional approach of solving the KS equation by an iterative 
algorithm, and that the dominant scaling in the calcula- 
tions of the energy functional and its derivative based on 
a plane-wave basis set is also 0(M 2 Nb as i S ) with Nf, as i S 
being the number of plane- wave basis functions, when ad- 
vantage is taken of the fast Fourier transform technique. 

In order to illustrate the utility of the DM method 
and the difficulty in the usage of the FD function, the 
electronic structures of C2, Si2 and a tungsten adatom 
adsorbed on the W(110) surface are calculated. Here- 
after, we obey the nine-point finite-difference formula 
for the derivative arising from the kinetic- energy opera- 
tor, and the dense-grid spacing is fixed at hdens = h/3 
H , where h denotes the coarse-grid spacing. The norm- 
conserving pseudopotential is employed in a Klcinman- 
Bylander nonlocal form pO 1 1 1 . 

We show in Figs. |l| and |^ the results of the application 
to the carbon dimer as a demonstration of the advantages 
of our method. It is well-known that a self-consistent 
screening potential for the carbon dimer cannot be sta- 
bilized by the integral occupation of the occupied states. 
Pederson and Jackson 0] showed by the study based on 
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the local-density approximation J12| that the minimum 
energy corresponds to a fractionally occupied system in 
the range of the atomic separation R—2A - 3.7 a.u. and 
that KS eigenvalues for the fractionally occupied states 
are degenerate. In Figs. |l] and ^, we took the grid spac- 
ing h=0. 33 a.u. in a cell of 16 3 a.u. under the nonperiodic 
boundary condition (b.c.) of vanishing wave function, 
and set the number of electrons N=8 and the number of 
wave functions M=5. Exchange-correlation effects were 
treated with the local-density approximation according 
to Pederson and Jackson Jl. Figure [l] illustrates the 
ground-state total energy, the KS eigenvalues of the lir u 
and 2a g states, and the occupation number of the 1tt u 
state as a function of atomic separation. The result of 
the brute-force method || implementing the direct vari- 
ation of the occupation numbers is also depicted here for 
comparison with our DM method. The total energies 
predicted by both methods coincide with one another 
[see Fig. 0(a)] . As shown in Fig. 0(b), KS eigenval- 
ues for the 1tt u and 2a g states are degenerate in the 
range i?=2.34 - 3.45 a.u. and the occupation numbers 
for these states become fractional in this range, which 
accords with the results of the other methods ||[||. In 
Fig. ^, we plot the occupation numbers of these states at 
R = 2.30 a.u. and 2.50 a.u. as a function of the num- 
ber of conjugate-gradient iterations, where the electronic 
structure at R = 2.40 a.u. is used as the starting point. 
One can see that the iterative process is markedly sta- 
ble and our DM procedure gives good convergence of the 
occupation numbers. 
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FIG. 1. C2 adiabatic-potential curve, eigenvalues of ln u 
and 2<r 9 states, lir u occupancy as a function of atomic sepa- 
ration R. In (b), the solid (dotted) curve represents the eigen- 
value of ln u (2<7 3 ) state, and the dash-dotted curve shows the 
occupation numebr of lir u state. Circles correspond to data 
obtained by the brute-force method. 
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Number of conjugate-gradient iterations 
FIG. 2. Occupation number n% = (2 — Af )Af as a function 
of the number of conjugate-gradient iterations. (For the defi- 
nition of A", see text preceding Eq. Q.) The solid (dotted) 
curves represent the occupation numbers in the case of the 
atomic separation R = 2.30 a.u. (2.50 a.u.); each run starts 
from the electronic structure at R = 2.40 a.u. 



We next give an example that the FD method includ- 
ing a fictitious broadening temperature leads to an incor- 
rect ground-state geometry in molecular-dynamics simu- 
lations. Figure |§| shows the curves of the thermodynamic 
free energy versus the atomic separation for the silicon 
dimer. The calculation was performed under the follow- 
ing conditions: the grid resolution h = 0.50 a.u., a cell of 
24 3 a.u. under the nonperiodic b.c. and the local-spin- 
density approximation Jl2]| for the exchange-correlation 
potential. The experiments (b|] proved that the energy 
difference between two molecular configurations (a) and 
(b) indicated in Fig. ^| is quite small, the equilibrium 
bond length is 4.07 a.u. for the configuration (a) and 4.25 
a.u. for (b), and the ground state with the lowest total 
energy is the latter. Some theoretical analyses have been 
carried out to examine the situation As seen 

in Fig. ||, our DM procedure can yield results that are 
in excellent agreement with the empirical data. On the 
contrary, the conventional KS scheme using the FD func- 
tion is unable to search out the minimum configuration 
(b) in molecular-dynamics iterations at the broadening 
temperature T>1 mH. 
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FIG. 3. Si2 adiabatic-potential curves as a function of 
atomic separation R obtained by the DM method and the con- 
ventional KS scheme at different broadening temperatures T. 
The indicated configurations are (a) (lcr g ) 2 (la u ) 2 (l7r u ) 4 and 
(b) (la 9 ) 2 (la 11 ) 2 (2a 9 ) 2 (l7r 11 ) 2 . 



3 



As a final example, we evaluate the energy barrier for a 
tungsten adatom to hop along the [111] direction on the 
W(110) surface. Our test system is a cell of 17 tungsten 
atoms with 102 electrons, consisting of one adatom and 
two rigid W(110) planes, under the periodic b.c. at the 
[001] and [110] directions, and the nonperiodic b.c. of 
vanishing wave function at the [110] one (i.e., thin film 
model) . The number of wave functions M — 64 for each 
spin, the coarse-grid spacing h = 0.30 a.u., and the local- 
spin-dcnsity approximation were employed. Only the V 
point of the Brillouin zone was sampled. We first eval- 
uate the ground-state geometry for the adatom located 
at point A in Fig. ^ using our DM approach, and then 
displace the adatom along the [111] direction from point 
A to B. The diffusion barrier of the tungsten adatom is 
given in Table |[ The numerical error in the barrier at 
T=2 mH is found to be not negligible but about 10 % of 
the true value obtained from our DM scheme. 



In conclusion, we have demonstrated a new variational 
technique to introduce appropriately the fractional oc- 
cupation numbers, and calculated the ground-state en- 
ergy and the charge density to a high degree of accuracy 
for a system with degeneracies at the Fermi level. Our 
scheme requires no additional statistical parameter such 
as a broadening temperature. Moreover, when a system 
is assumed to be described in terms of localized wave 
functions, our procedure makes it possible to perform 
the calculation with linear system-size scaling O(N). Re- 
search in this direction is in progress. 

One of the authors (T.O.) is a Research Fellow of the 
Japan Society for the Promotion of Science. This work 
was partially supported by a Grant-in- Aid for COE Re- 
search. The numerical calculation was carried out by the 
computer facilities at the Institute for Solid State Physics 
at the University of Tokyo. 
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FIG. 4. Top view of W(110) surface-layer atoms (large 
open circles), second-layer atoms (small open circles) and 
adatom (crosses) for the surface diffusion discussed in the 
text. 



TABLE I. The diffusion barrier of a tungsten adatom on 
the rigid W(110) surface along the [111] direction. Data from 
our DM method and the KS scheme including a broadening 
temperature are presented. 



Temperature (Hartree) 


Diffusion barrier (eV) 


Present work 


1.40 


0.001 


1.43 


0.002 


1.54 


0.004 


1.54 
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